function dif_file( idlPath, idlFile, strPath, strFile )
%--------------------------------------------------------------------------
% Error Spectrogram of Audio File
%
% Cooper Baker - 2014
%--------------------------------------------------------------------------

close all;

% Settings
%--------------------------------------------------------------------------
winSize = 512;
overlap = 4;
window  = chebwin( winSize, 125 );
name    = 'Error Spectrogram';

% Initializations
%--------------------------------------------------------------------------
if any( exist( 'idlFile' ) ~= 1 )
    [ idlFile, idlPath ] = uigetfile( '*.wav', 'Ideal Audio File' );
    [ strFile, strPath ] = uigetfile( '*.wav', 'Stretched Audio File' );
end

[ idlBuf, sr ] = audioread( [ idlPath, idlFile ] );
[ strBuf, sr ] = audioread( [ strPath, strFile ] );
idlSize        = length( idlBuf );
strSize        = length( strBuf );
halfWinSize    = winSize / 2;
overlapSamps   = winSize - ( winSize / overlap );
hopSize        = winSize / overlap;
hopMax         = length( idlBuf ) - winSize;
idlMsec        = ( idlSize / sr ) * 1000;
gridMsec       = 500;
gridHops       = gridMsec / ( ( ( winSize / overlap ) / sr ) * 1000 );
hopMsec        = ( hopSize / sr ) * 1000;
index          = 1;
winSum         = sum( window );

% crop or pad stretch buffer to match ideal buffer times two
if ( idlSize * 2 ) > strSize
    strBuf = [ strBuf ; zeros( ( idlSize * 2 ) - strSize, 1 ) ];
elseif strSize > ( idlSize * 2 )
    strBuf = strBuf( 1 : idlSize * 2 );
end

strSize = length( strBuf );

% Normalization
%--------------------------------------------------------------------------

% get file metadata
idlInfo = audioinfo( [ idlPath, idlFile ] );
strInfo = audioinfo( [ strPath, strFile ] );

% get file normalization coefficient
idlFileNorm = str2num( idlInfo.Title );
strFileNorm = str2num( strInfo.Title );

% de-normalize normalized .wav file data
idlRaw = idlBuf / idlFileNorm;
strRaw = strBuf / strFileNorm;

% make hann windows
idlWin = hann( idlSize );
strWin = hann( strSize );

% copy and window files
idlNormWin = idlRaw .* idlWin;
strNormWin = strRaw .* strWin;

% calculate rms of windowed file middles
idlRms = rms( idlNormWin );
strRms = rms( strNormWin );

% calculate normalization scalar
strScale = 1 / ( strRms / idlRms );

% scale stretch buffer to match ideal buffer
idlNorm = idlRaw;
strNorm = strRaw * strScale;

% Analysis
%--------------------------------------------------------------------------
idlGram     = spectrogram( idlNorm, window, overlapSamps, winSize, sr );
strGramLong = spectrogram( strNorm, window, overlapSamps, winSize, sr );

idlSize = length( idlGram ) + 1;

% compress stretched array to match ideal array
while index < idlSize;

    % calculate stretched array index
    strIndex = ( index * 2 - 1 );

    % average adjacent pairs of spectra
    strGram( :, index ) = ( strGramLong( :, strIndex ) + strGramLong( :, strIndex + 1 ) ) ./ 2;

    %increment index
    index = index + 1;
end

% calculate magnitude spectra
idlMag = abs( idlGram );
strMag = abs( strGram );

% calculate amplitude spectra
idlAmp = ( idlMag / winSum );
strAmp = ( strMag / winSum );

% calculate sones spectra
idlSones = idlAmp .^ 0.6;
strSones = strAmp .^ 0.6;

% copy spectra to plot arrays
idlPlot = idlSones;
strPlot = strSones;

% calculate ideal plot offset
offset = 0;
offset = abs( min( [ min( idlPlot ) min( strPlot ) ] ) );

% offset plots
idlPlot = idlPlot + offset;
strPlot = strPlot + offset;

% calculate absolute difference between plots (error)
diffPlot = abs( idlPlot - strPlot );

% find maximum value of difference plot
diffMax = max( diffPlot(:) );

% Plot
%--------------------------------------------------------------------------
fontName   = 'Times New Roman';
fontSize   = 12;

fig = figure( 1 );
set( fig, 'Name', name );
set( fig, 'Position', [ 0 0 800 250 ] );
set( fig, 'defaultAxesFontName', fontName );
set( fig, 'defaultTextFontName', fontName );

colormap( [ [ 1 : -0.2/64 : 0.8 ]; ( [ 1 : -0.5/64 : 0.5 ] .^ 10 ) ; [ 1 : -0.5/64 : 0.5 ] .^ 20 ]' );

% spectrogram
imagesc( diffPlot );
caxis( 'auto' );
axis xy;

% grid and labels
axis( [ 0 length( diffPlot ) 1 halfWinSize ] );
set( gca, 'XTick', [ 0 : gridHops : length( diffPlot ) ] );
set( gca, 'XTickLabel', [ 0 : gridMsec : idlMsec ] );
set( gca, 'YTick', [ 1 : ( halfWinSize / ( sr / 4000 ) ) : halfWinSize ] );
set( gca, 'YTickLabel', [ 0 : 2000 : ( sr / 2 ) ] );
xlabel( 'Milliseconds' );
ylabel( 'Frequency' );

% legend
bar = colorbar( 'location', 'East' );
set( get( bar, 'YLabel'), 'String', 'Sones' );
set( get( bar, 'YLabel'), 'Position', [ 0.75 ( diffMax * 0.1 ) ] );

% tighten up figure borders
tightfig();

% Save Data to Files
%--------------------------------------------------------------------------

% get file name
[ a, fileName, b ] = fileparts( strFile );

% write plot to file
hgexport( fig, [ fileName, '.spect.eps' ] );

% write error data to files
csvFile = sprintf( '%s.spect.csv', fileName );
csvwrite( csvFile, diffPlot );

% EOF
